Properties of the Fourier Transform¶

These demonstrations go beyond the content of the course. I don't necessarily expect you to understand all that is going on in these snippets of code. But I show them to you so you see what sorts of things can be implemented by manipulating the Fourier coefficients.

Preliminaries¶

In [1]:
import numpy as np
from numpy.fft import fft, ifft, fft2, ifft2, fftshift, ifftshift
import matplotlib.pyplot as plt
from scipy import ndimage, misc
from scipy.signal import gaussian
In [2]:
from IPython.display import display, HTML
display(HTML("<style>.container { width:98% !important; }</style>"))

np.set_printoptions(edgeitems=30, linewidth=250, 
    formatter=dict(float=lambda x: "% .4f" % x))
In [3]:
def Gaussian2D(dims, std):
    '''
        Creates an array containing a 2-D Gaussian function, with mean
        at mu, and standard deviation std.

        Inputs:
          dims is a 2-vector containing the number of rows and columns
          std is the standard deviation
    '''
    nr, nc = dims
    rr, cc = np.mgrid[0:nr,0:nc]
    ctr = GetCentre(rr)
    temp1 = rr - ctr[0]
    temp2 = cc - ctr[1]
    temp1 = np.exp(-temp1**2/std)
    temp2 = np.exp(-temp2**2/std)
    blah = temp1*temp2
    return blah/sum(blah.flatten())

def GetCentre(f):
    return np.array(np.floor(np.array(np.shape(f))/2.), dtype=int)

Conjugate Symmetry¶

In [4]:
# If f is real-valued...
f = np.random.random(8)
print(f)

[ 0.8983  0.3727  0.9525  0.8016  0.4681  0.1283  0.3566  0.1619]
In [5]:
# Fourier coefs will be cojugate symmetric
F = np.fft.fft(f)
print(F)
print(F[2], F[-2])

[4.13988582+0.j         0.15064761-1.22107551j 0.05739223+0.46241617j 0.70967138-0.0292481j  1.21095411+0.j         0.70967138+0.0292481j  0.05739223-0.46241617j 0.15064761+1.22107551j]
(0.057392232805012044+0.46241617323556405j) (0.057392232805012044-0.46241617323556405j)

Convolution¶

In [6]:
rr, cc = np.mgrid[-128:128, -128:128]

Try a circular kernel¶

In [7]:
d = np.sqrt(rr**2 + cc**2)
circle = np.zeros([256,256])
circle_mask = d<15
circle[circle_mask] = 1.
circle = ifftshift(circle)
circle[0:22,0:5] = 1.
plt.figure(figsize=(8,4))
plt.subplot(1,2,1); plt.imshow(fftshift(circle), cmap='gray'); plt.axis('on'); plt.title('g, Shifted');
plt.subplot(1,2,2); plt.imshow(circle, cmap='gray'); plt.axis('on'); plt.title('g, Unshifted');
In [8]:
# Create an image with just a few non-zero pixels.
f = np.zeros([256,256])
f[50,100] = 1.0
f[48,110] = 2.
plt.figure(figsize=(8,8))
plt.imshow(f, cmap='gray'); plt.axis('off'); plt.title('f');
In [9]:
g = ifft2( fft2(f) * fft2(circle))

plt.figure(figsize=(8,8))
plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title(r'$(f \ast g)$');
In [10]:
f = plt.imread('pd.jpg')
f = f[:,:,0]
plt.imshow(f, cmap='gray'); plt.axis('off'); plt.title('An image f');
In [11]:
g = ifft2( fft2(f) * fft2(circle))

#plt.figure(figsize=(8,8))
plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title('g convolved with the image');

Or, try a Gaussian edge-detector kernel¶

In [12]:
# Make and edge detector
g = Gaussian2D([256,256], 20)
G = ifftshift( fft2(fftshift(g)) )
ramp = np.outer(range(256), range(256)) - 128**2
H = G * ramp*1.j
circle = np.real( ifft2( fftshift(H) ) )
plt.figure(figsize=[8,5]);
plt.subplot(1,2,1); plt.imshow(fftshift(circle), cmap='gray'); plt.axis('off'); plt.title('Edge Kernel (Shifted)');
plt.subplot(1,2,2); plt.imshow(circle, cmap='gray'); plt.axis('off'); plt.title('Edge Kernel (Unshifted)');
In [13]:
f = plt.imread('pd.jpg')
f = f[:,:,0]
In [14]:
F = fftshift( fft2( f ) )
In [15]:
plt.figure(figsize=[10,20])
plt.subplot(1,2,1)
plt.imshow(np.log(abs(F)+1), cmap='gray'); plt.axis('off'); plt.title('Shifted DFT (log-modulus)');
plt.subplot(1,2,2)
plt.imshow(f, cmap='gray'); plt.axis('off');
In [16]:
g = fftshift( ifft2( fft2(ifftshift(f)) * fft2(circle)) )

plt.figure(figsize=(16,8))
plt.subplot(1,2,1);
G = fftshift(fft2(ifftshift(g)))
plt.imshow(np.abs(G), cmap='gray'); plt.axis('off')
plt.title('DFT(Edge Kernel) x DFT(Image)');
plt.subplot(1,2,2); plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title('Edge Kernel convolved with image');

Filtering in the frequency domain¶

In [17]:
thresh = 10
rows, cols = np.shape(f)
ctr = np.floor(np.array(np.shape(f))/2) +1
rr, cc = np.mgrid[-ctr[0]:(rows-ctr[0]),
                  -ctr[1]:(cols-ctr[1])]
In [18]:
def PlotFiltImg(F, filt):
    F_filt = (filt)*F
    f_filt = ifft2( fftshift(F_filt) )
    plt.figure(figsize=[15,15])
    plt.subplot(2,2,1); plt.imshow(f, cmap='gray'); plt.axis('off');
    plt.subplot(2,2,3); plt.imshow(np.log(abs(F)+1), cmap='gray'); plt.axis('off');
    plt.subplot(2,2,4); plt.imshow(np.log(abs(F_filt)+1), cmap='gray'); plt.axis('off');
    plt.subplot(2,2,2); plt.imshow(np.real(f_filt), cmap='gray'); plt.axis('off');
In [19]:
filt_mask = abs(cc)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
In [20]:
filt_mask = abs(rr)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
In [21]:
filt_mask = np.sqrt(rr**2 + cc**2)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
In [22]:
# Make and edge detector
g = Gaussian2D([256,256], 5)
G = ifftshift( fft2(fftshift(g)) )
# ramp = np.outer(range(256), np.ones(256)) - 128**2
# Hv = G * ramp*1.j
# ramp = np.outer(np.ones(256), range(256)) - 128**2
# Hh = G * ramp*1.j

Applying that filter in the frequency domain is the same as convolving the image with the kernel below, attained by taking the IDFT of the filter kernel.

In [23]:
Filt = fftshift(ifft2(ifftshift(filt)))
plt.figure()
plt.subplot(1,2,1); plt.imshow(filt, cmap='gray'); plt.title('DFT of filter'); plt.axis('off')
plt.subplot(1,2,2); plt.imshow(np.log(abs(Filt)+1), cmap='gray'); plt.axis('off'); plt.title('Spatial filter');

Scaling¶

You can scale an image using the DFT by padding in cropping in the two domains (time/space, and frequency).

In [24]:
f = plt.imread('pd.jpg')
f = f[:,:,0]

# We will magnify an image by padding its DFT coefficients, and then cropping the image
# that results from the IDFT of those padded coefs.
scale_factor = 1.5
padamt = int(rows*(scale_factor-1.)/2.)
G = np.pad(fftshift(fft2(f)), (padamt,), mode='constant')
g512 = ifft2( ifftshift(G) )
if padamt!=0:
    g = g512[padamt:-padamt,padamt:-padamt]
else:
    g = g512
In [25]:
plt.figure(1,figsize=[15,15])
plt.clf()
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray'); plt.title('Scaled x'+str(scale_factor));

Fourier Shift Theorem¶

You can shift an image by adding a linear 'ramp' to the phase portion of its DFT coefficients.

In [38]:
dr = 4
dc = 12
ramp = -2.*np.pi*1j*(cc*dc + rr*dr)/rows
ramp = np.roll(np.roll(ramp,-1, axis=0),-1, axis=1)
wave = np.exp(ramp)
plt.figure(figsize=(10,5))
plt.subplot(1,2,1); plt.imshow(np.real(wave))
plt.subplot(1,2,2); plt.imshow(np.imag(ramp))
plt.title('Ramp in imaginary part');
In [39]:
# Make sure the centre pixel of the ramp has a phase of zero.
ctr = int(np.floor(256/2))
print('This should be zero -> '+str(ramp[ctr,ctr]))
print('This is the DC of our image -> '+str(F[ctr,ctr]))  # This should be the DC of our image.

This should be zero -> -0j
This is the DC of our image -> (2491737+0j)
In [40]:
G = wave * fftshift(fft2(f))
g = ifft2(ifftshift(G))

plt.figure(1,figsize=[15,15])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray');
plt.title('Shifted (down,right)=('+str(dr)+','+str(dc)+')');

Pendulum synchronization video (YouTube)

Phase Correlation¶

You can detect a shift between two similar images by looking for a ramp in their phase difference. Or, you can take that phase difference and take the IDFT of it, and look for a spike. The offset of the spike from the origin gives the relative shift.

In [54]:
plt.figure(figsize=(8,8));
g = ndimage.shift(f,[dr,dc]) + np.random.normal(0.0, 5, size=np.shape(f))
G = fft2(g)
F = fft2(f)
H = G/F
plt.imshow(np.real(fftshift(H)), cmap='gray', vmin=-4, vmax=4);
In [55]:
plt.figure(figsize=(8,8));
h = ifft2(H)
plt.imshow(np.real(h), cmap='gray');
In [56]:
# Where does that lonely pixel occur?
coords = np.unravel_index(np.argmax(np.real(h)), h.shape)
print('Bright pixel at '+str(coords))
print('True shift was '+str((dr,dc)))

Bright pixel at (4, 12)
True shift was (4, 12)

Rotation¶

You can also rotate an image by simply applying the same rotation to the frequency domain.

In [41]:
theta = 30
f = plt.imread('pd.jpg')
f = np.array(f[:,:,0], dtype=float)
In [47]:
padamt = int(rows/2)
#padamt = 0
In [48]:
# Padding the image first makes the frequency domain smoother
ff = f
ff = np.pad(f, (padamt,), mode='constant')
FF = fftshift( fft2( ifftshift(ff) ) )
In [49]:
# Apply rotation to the Fourier coefficients
Gr = ndimage.rotate(np.real(FF), theta, reshape=False, order=1)
Gi = ndimage.rotate(np.imag(FF), theta, reshape=False, order=1)
G = Gr + 1j*Gi
In [50]:
g = fftshift( ifft2( ifftshift(G) ) )
if padamt!=0:
    g = g[padamt:-padamt,padamt:-padamt]
In [51]:
plt.figure(figsize=[15,15])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original'); plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray'); 
plt.title('Rotated by '+str(theta)+' degrees'); plt.axis('off');

Fourier Projection Theorem¶

The Fourier projection theorem is the basis for reconstruction in Computed Tomography (CT) scans.

In [57]:
f = plt.imread('ct.jpg')
f1 = np.array(f[:,:,0], dtype=float)
In [58]:
plt.figure(1, figsize=(10,10))
plt.clf()
plt.subplot(2,2,1)
plt.imshow(f1, cmap=plt.cm.gray);

px = np.sum(f1,axis=0)
plt.subplot(2,2,3)
plt.plot(px)
plt.axis('tight');

F = fftshift(fft2(f1))
plt.subplot(2,2,2)
Fnice = np.log(abs(F)+1)
plt.imshow(Fnice, cmap=plt.cm.gray);

ctr = GetCentre(F)
Fx = F[ctr[0],:]
plt.plot([1,254],[128,128], 'r--')
plt.axis([0,255,255,0])

plt.subplot(2,2,4)
plt.plot(np.log(abs(Fx)+1), 'r')
plt.axis('tight');

Px = fftshift( np.fft.fft(px) )
#plt.subplot(2,2,4)
plt.plot(np.log(abs(Px)+1), 'b--');
In [ ]: